home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / sinint.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  10.6 KB  |  403 lines
  1. /* specfunc/sinint.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_trig.h>
  26. #include <gsl/gsl_sf_expint.h>
  27.  
  28. #include "error.h"
  29.  
  30. #include "chebyshev.h"
  31. #include "cheb_eval.c"
  32.  
  33. /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  34.  
  35. /* based on SLATEC r9sifg.f, W. Fullerton */
  36.  
  37. /*
  38.  series for f1   on the interval  2.00000e-02 to  6.25000e-02
  39.                     with weighted error   2.82e-17
  40.                      log weighted error  16.55
  41.                    significant figures required  15.36
  42.                     decimal places required  17.20
  43. */
  44. static double f1_data[20] = {
  45.    -0.1191081969051363610,
  46.    -0.0247823144996236248,
  47.     0.0011910281453357821,
  48.    -0.0000927027714388562,
  49.     0.0000093373141568271,
  50.    -0.0000011058287820557,
  51.     0.0000001464772071460,
  52.    -0.0000000210694496288,
  53.     0.0000000032293492367,
  54.    -0.0000000005206529618,
  55.     0.0000000000874878885,
  56.    -0.0000000000152176187,
  57.     0.0000000000027257192,
  58.    -0.0000000000005007053,
  59.     0.0000000000000940241,
  60.    -0.0000000000000180014,
  61.     0.0000000000000035063,
  62.    -0.0000000000000006935,
  63.     0.0000000000000001391,
  64.    -0.0000000000000000282
  65. };
  66. static cheb_series f1_cs = {
  67.   f1_data,
  68.   19,
  69.   -1, 1,
  70.   10
  71. };
  72.  
  73. /*
  74.  
  75.  series for f2   on the interval  0.00000e+00 to  2.00000e-02
  76.                     with weighted error   4.32e-17
  77.                      log weighted error  16.36
  78.                    significant figures required  14.75
  79.                     decimal places required  17.10
  80. */
  81. static double f2_data[29] = {
  82.    -0.0348409253897013234,
  83.    -0.0166842205677959686,
  84.     0.0006752901241237738,
  85.    -0.0000535066622544701,
  86.     0.0000062693421779007,
  87.    -0.0000009526638801991,
  88.     0.0000001745629224251,
  89.    -0.0000000368795403065,
  90.     0.0000000087202677705,
  91.    -0.0000000022601970392,
  92.     0.0000000006324624977,
  93.    -0.0000000001888911889,
  94.     0.0000000000596774674,
  95.    -0.0000000000198044313,
  96.     0.0000000000068641396,
  97.    -0.0000000000024731020,
  98.     0.0000000000009226360,
  99.    -0.0000000000003552364,
  100.     0.0000000000001407606,
  101.    -0.0000000000000572623,
  102.     0.0000000000000238654,
  103.    -0.0000000000000101714,
  104.     0.0000000000000044259,
  105.    -0.0000000000000019634,
  106.     0.0000000000000008868,
  107.    -0.0000000000000004074,
  108.     0.0000000000000001901,
  109.    -0.0000000000000000900,
  110.     0.0000000000000000432
  111. };
  112. static cheb_series f2_cs = {
  113.   f2_data,
  114.   28,
  115.   -1, 1,
  116.   14
  117. };
  118.  
  119. /*
  120.  
  121.  series for g1   on the interval  2.00000e-02 to  6.25000e-02
  122.                     with weighted error   5.48e-17
  123.                      log weighted error  16.26
  124.                    significant figures required  15.47
  125.                     decimal places required  16.92
  126. */
  127. static double g1_data[21] = {
  128.    -0.3040578798253495954,
  129.    -0.0566890984597120588,
  130.     0.0039046158173275644,
  131.    -0.0003746075959202261,
  132.     0.0000435431556559844,
  133.    -0.0000057417294453025,
  134.     0.0000008282552104503,
  135.    -0.0000001278245892595,
  136.     0.0000000207978352949,
  137.    -0.0000000035313205922,
  138.     0.0000000006210824236,
  139.    -0.0000000001125215474,
  140.     0.0000000000209088918,
  141.    -0.0000000000039715832,
  142.     0.0000000000007690431,
  143.    -0.0000000000001514697,
  144.     0.0000000000000302892,
  145.    -0.0000000000000061400,
  146.     0.0000000000000012601,
  147.    -0.0000000000000002615,
  148.     0.0000000000000000548
  149. };
  150. static cheb_series g1_cs = {
  151.   g1_data,
  152.   20,
  153.   -1, 1,
  154.   13
  155. };
  156.  
  157. /*
  158.  
  159.  series for g2   on the interval  0.00000e+00 to  2.00000e-02
  160.                     with weighted error   5.01e-17
  161.                      log weighted error  16.30
  162.                    significant figures required  15.12
  163.                     decimal places required  17.07
  164. */
  165. static double g2_data[34] = {
  166.    -0.0967329367532432218,
  167.    -0.0452077907957459871,
  168.     0.0028190005352706523,
  169.    -0.0002899167740759160,
  170.     0.0000407444664601121,
  171.    -0.0000071056382192354,
  172.     0.0000014534723163019,
  173.    -0.0000003364116512503,
  174.     0.0000000859774367886,
  175.    -0.0000000238437656302,
  176.     0.0000000070831906340,
  177.    -0.0000000022318068154,
  178.     0.0000000007401087359,
  179.    -0.0000000002567171162,
  180.     0.0000000000926707021,
  181.    -0.0000000000346693311,
  182.     0.0000000000133950573,
  183.    -0.0000000000053290754,
  184.     0.0000000000021775312,
  185.    -0.0000000000009118621,
  186.     0.0000000000003905864,
  187.    -0.0000000000001708459,
  188.     0.0000000000000762015,
  189.    -0.0000000000000346151,
  190.     0.0000000000000159996,
  191.    -0.0000000000000075213,
  192.     0.0000000000000035970,
  193.    -0.0000000000000017530,
  194.     0.0000000000000008738,
  195.    -0.0000000000000004487,
  196.     0.0000000000000002397,
  197.    -0.0000000000000001347,
  198.     0.0000000000000000801,
  199.    -0.0000000000000000501
  200. };
  201. static cheb_series g2_cs = {
  202.   g2_data,
  203.   33,
  204.   -1, 1,
  205.   20
  206. };
  207.  
  208.  
  209. /* x >= 4.0 */
  210. static void fg_asymp(const double x, gsl_sf_result * f, gsl_sf_result * g)
  211. {
  212.   /*
  213.       xbig = sqrt (1.0/r1mach(3))
  214.       xmaxf = exp (amin1(-alog(r1mach(1)), alog(r1mach(2))) - 0.01)
  215.       xmaxg = 1.0/sqrt(r1mach(1))
  216.       xbnd = sqrt(50.0)
  217.   */
  218.   const double xbig  = 1.0/GSL_SQRT_DBL_EPSILON;
  219.   const double xmaxf = 1.0/GSL_DBL_MIN;
  220.   const double xmaxg = 1.0/GSL_SQRT_DBL_MIN;
  221.   const double xbnd  = 7.07106781187;
  222.  
  223.   const double x2 = x*x;
  224.  
  225.   if(x <= xbnd) {
  226.     gsl_sf_result result_c1;
  227.     gsl_sf_result result_c2;
  228.     cheb_eval_e(&f1_cs, (1.0/x2-0.04125)/0.02125, &result_c1);
  229.     cheb_eval_e(&g1_cs, (1.0/x2-0.04125)/0.02125, &result_c2);
  230.     f->val = (1.0 + result_c1.val)/x;
  231.     g->val = (1.0 + result_c2.val)/x2;
  232.     f->err = result_c1.err/x  + 2.0 * GSL_DBL_EPSILON * fabs(f->val);
  233.     g->err = result_c2.err/x2 + 2.0 * GSL_DBL_EPSILON * fabs(g->val);
  234.   }
  235.   else if(x <= xbig) {
  236.     gsl_sf_result result_c1;
  237.     gsl_sf_result result_c2;
  238.     cheb_eval_e(&f2_cs, 100.0/x2-1.0, &result_c1);
  239.     cheb_eval_e(&g2_cs, 100.0/x2-1.0, &result_c2);
  240.     f->val = (1.0 + result_c1.val)/x;
  241.     g->val = (1.0 + result_c2.val)/x2;
  242.     f->err = result_c1.err/x  + 2.0 * GSL_DBL_EPSILON * fabs(f->val);
  243.     g->err = result_c2.err/x2 + 2.0 * GSL_DBL_EPSILON * fabs(g->val);
  244.   }
  245.   else {
  246.     f->val = (x < xmaxf ? 1.0/x  : 0.0);
  247.     g->val = (x < xmaxg ? 1.0/x2 : 0.0);
  248.     f->err = 2.0 * GSL_DBL_EPSILON * fabs(f->val);
  249.     g->err = 2.0 * GSL_DBL_EPSILON * fabs(g->val);
  250.   }
  251.  
  252.   return;
  253. }
  254.  
  255.  
  256. /* based on SLATEC si.f, W. Fullerton
  257.  
  258.  series for si   on the interval  0.00000e+00 to  1.60000e+01
  259.                     with weighted error   1.22e-17
  260.                      log weighted error  16.91
  261.                    significant figures required  16.37
  262.                     decimal places required  17.45
  263. */
  264.  
  265. static double si_data[12] = {
  266.   -0.1315646598184841929,
  267.   -0.2776578526973601892,
  268.    0.0354414054866659180,
  269.   -0.0025631631447933978,
  270.    0.0001162365390497009,
  271.   -0.0000035904327241606,
  272.    0.0000000802342123706,
  273.   -0.0000000013562997693,
  274.    0.0000000000179440722,
  275.   -0.0000000000001908387,
  276.    0.0000000000000016670,
  277.   -0.0000000000000000122
  278. };
  279.  
  280. static cheb_series si_cs = {
  281.   si_data,
  282.   11,
  283.   -1, 1,
  284.   9
  285. };
  286.  
  287. /*
  288.  series for ci   on the interval  0.00000e+00 to  1.60000e+01
  289.                     with weighted error   1.94e-18
  290.                      log weighted error  17.71
  291.                    significant figures required  17.74
  292.                     decimal places required  18.27
  293. */
  294. static double ci_data[13] = {
  295.    -0.34004281856055363156,
  296.    -1.03302166401177456807,
  297.     0.19388222659917082877,
  298.    -0.01918260436019865894,
  299.     0.00110789252584784967,
  300.    -0.00004157234558247209,
  301.     0.00000109278524300229,
  302.    -0.00000002123285954183,
  303.     0.00000000031733482164,
  304.    -0.00000000000376141548,
  305.     0.00000000000003622653,
  306.    -0.00000000000000028912,
  307.     0.00000000000000000194
  308. };
  309. static cheb_series ci_cs = {
  310.   ci_data,
  311.   12,
  312.   -1, 1,
  313.   9
  314. };
  315.  
  316.  
  317. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  318.  
  319. int gsl_sf_Si_e(const double x, gsl_sf_result * result)
  320. {
  321.   double ax = fabs(x);
  322.   
  323.   /* CHECK_POINTER(result) */
  324.  
  325.   if(ax < GSL_SQRT_DBL_EPSILON) {
  326.     result->val = x;
  327.     result->err = 0.0;
  328.     return GSL_SUCCESS;
  329.   }
  330.   else if(ax <= 4.0) {
  331.     gsl_sf_result result_c;
  332.     cheb_eval_e(&si_cs, (x*x-8.0)*0.125, &result_c);
  333.     result->val  =  x * (0.75 + result_c.val);
  334.     result->err  = ax * result_c.err;
  335.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  336.     return GSL_SUCCESS;
  337.   }
  338.   else {
  339.     /* Note there is no loss of precision
  340.      * here bcause of the leading constant.
  341.      */
  342.     gsl_sf_result f;
  343.     gsl_sf_result g;
  344.     fg_asymp(ax, &f, &g);
  345.     result->val  = 0.5 * M_PI - f.val*cos(ax) - g.val*sin(ax);
  346.     result->err  = f.err + g.err;
  347.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  348.     if(x < 0.0) result->val = -result->val;
  349.     return GSL_SUCCESS;
  350.   }
  351. }
  352.  
  353.  
  354. int gsl_sf_Ci_e(const double x, gsl_sf_result * result)
  355. {
  356.   /* CHECK_POINTER(result) */
  357.  
  358.   if(x <= 0.0) {
  359.     DOMAIN_ERROR(result);
  360.   }
  361.   else if(x <= 4.0) {
  362.     const double lx = log(x);
  363.     const double y  = (x*x-8.0)*0.125;
  364.     gsl_sf_result result_c;
  365.     cheb_eval_e(&ci_cs, y, &result_c);
  366.     result->val  = lx - 0.5 + result_c.val;
  367.     result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(lx) + 0.5) + result_c.err;
  368.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  369.     return GSL_SUCCESS;
  370.   }
  371.   else {
  372.     gsl_sf_result sin_result;
  373.     gsl_sf_result cos_result;
  374.     int stat_sin = gsl_sf_sin_e(x, &sin_result);
  375.     int stat_cos = gsl_sf_cos_e(x, &cos_result);
  376.     gsl_sf_result f;
  377.     gsl_sf_result g;
  378.     fg_asymp(x, &f, &g);
  379.     result->val  = f.val*sin_result.val - g.val*cos_result.val;
  380.     result->err  = fabs(f.err*sin_result.val);
  381.     result->err += fabs(g.err*cos_result.val);
  382.     result->err += fabs(f.val*sin_result.err);
  383.     result->err += fabs(g.val*cos_result.err);
  384.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  385.     return GSL_ERROR_SELECT_2(stat_sin, stat_cos);
  386.   }
  387. }
  388.  
  389.  
  390. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  391.  
  392. #include "eval.h"
  393.  
  394. double gsl_sf_Si(const double x)
  395. {
  396.   EVAL_RESULT(gsl_sf_Si_e(x, &result));
  397. }
  398.  
  399. double gsl_sf_Ci(const double x)
  400. {
  401.   EVAL_RESULT(gsl_sf_Ci_e(x, &result));
  402. }
  403.